
function [ceqNum , ceq , Gceq] = Bcr4bp_StateConsEMRot2D_Coe(tau_k , xx_tauk , var_tauk , ...
    tau_kp1 , xx_taukp1 , var_taukp1 , ...
    xx_taukp1_ , phi_tauk_taukp1 , nVar , nNode , jLoop , aux)
%
% 构造状态连接约束
%
% 注意：
% 如果xxkp1是'dv'，则约束为r1f - r20
% 否则都为x1f - x20
%
% 作者：张晨
% 邮箱：chenzhang@csu.ac.cn
% 单位：中国科学院空间应用工程与技术中心，空间探索室
% 时间：2021年08月25日
%%%%%%%%%%%%%%%%%%%%%%%%

% ---------------------------------------------------------------------------
% tauk动力学
dxx_tauk = Bcr4bp_EqmEMRot2D(tau_k , xx_tauk , aux);

% taukp1_动力学
dxx_taukp1_ = Bcr4bp_EqmEMRot2D(tau_kp1 , xx_taukp1_ , aux);

% 模块
dtauj_dtau1 = 1 - (jLoop - 1) / (nNode - 1);
dtaujp1_dtau1 = 1 - jLoop / (nNode - 1);
dtauj_dtaun = (jLoop - 1) / (nNode - 1);
dtaujp1_dtaun = jLoop / (nNode - 1);

% ---------------------------------------------------------------------------
% 连接约束数
ceqNum = 4;

% 约束
ceq = (xx_taukp1_ - xx_taukp1);

% 梯度预分配
Gceq = zeros(4 , nVar);

% --------------- 状态连接约束（梯度） -------------
% wrt xx_tauk
if jLoop == 1

    Gceq(: , 4 * (jLoop - 1) + 1 : 4 * jLoop) = phi_tauk_taukp1  * dxde_leo(var_tauk , aux);

else

    Gceq(: , 4 * (jLoop - 1) + 1 : 4 * jLoop) = phi_tauk_taukp1;

end

% wrt xxkp1
if jLoop == (nNode - 1)

    Gceq(: , 4 * jLoop + 1 : 4 * (jLoop + 1)) = - eye(4) * dxde_llo(var_taukp1 , aux);

else

    Gceq(: , 4 * jLoop + 1 : 4 * (jLoop + 1)) = - eye(4);

end

% 状态连接约束 wrt tau_1
Gceq(: , 4 * nNode + 1) = - phi_tauk_taukp1 * dxx_tauk * dtauj_dtau1 + dxx_taukp1_ * dtaujp1_dtau1;

% 状态连接约束 wrt tau_n
Gceq(: , 4 * nNode + 2) = - phi_tauk_taukp1 * dxx_tauk * dtauj_dtaun + dxx_taukp1_ * dtaujp1_dtaun;

end
